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ABSTRACT 

Aims. The large quantity and high quality of modern radio and infrared line observations require efficient modeling techniques 
to infer physical and chemical parameters such as temperature, density, and molecular abundances. 

Methods. We present a computer program to calculate the intensities of atomic and molecular lines produced in a uniform 
medium, based on statistical equilibrium calculations involving collisional and radiative processes and including radiation from 
background sources. Optical depth effects are treated with an escape probability method. The program is available on the World 
Wide Web at http : / /www. sron. rug . nl/~vdtak/radex/index. shtml . The program makes use of molecular data files 
maintained in the Leiden Atomic and Molecular Database (LAMDA), which will continue to be improved and expanded. 
Results. The performance of the program is compared with more approximate and with more sophisticated methods. An 
Appendix provides diagnostic plots to estimate physical parameters from line intensity ratios of commonly observed molecules. 
Conclusions. This program should form an important tool in analyzing observations from current and future radio and infrared 
telescopes. 
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1. Introduction 

Observations of spectral lines at radio, (sub)millimeter 
and infrared wavelengths are a powerful tool to in- 
vestigate the physical and chemical conditions in the 
dilute gas of astronomical sources where thermody- 



namic equilibrium is a poor approximation (e.g. JGenzel 



Il99ll Black 2000). To extract astrophysical parame- 
ters from the data, the excitation and optical depth of 
the lines need to be estimated, for which various meth- 
ods m ay be used, depending on the available observa- 
tions dVan Dishoeck & Hogerheijddll999l: IVan der Tald 



20051) . 

If only one or two lines of a moleculd^ have been ob- 
served, the excitation must be deduced from observa- 
tions of other species or from theoretical considerations. 
An example is the assumption that the excitation tem- 
perature equals the kinetic temperature, a case known as 
Local Thermodynamic Equilibrium (LTE) which holds 
at high densities. 



If many lines have been observed, a popular method 
is the 'rotation diagram', also c alled 'Boltzmann 
plot' or 'population diagram' (e.g.. [Blake et al.lll987l: 
Helmich et ail 1 1994 iGoldsmifh & Langerlll999f) . This 
method describes the excitation by a single tempera- 
ture, obtained by a fit to the line intensities as a func- 
tion of upper level energy. Provided that beam sizes are 
similar and optical depths are low, or that appropriate 
corrections are made, this method yields estimates of 
the excitation temperature and column density of the 
molecule. The excitation temperature approaches the ki- 
netic temperature in the high-density limit, but generally 
depends on both kinetic temperature and volume den- 
sity. Spectral line surveys are often analyzed with ro- 
tation diag rams, although more advanced methods are 
also u sed ( Helmich & van Dishoeck! 1997 ; Comito et al 
2005). 



More sophisticated methods retain the assumption of 
a local excitation, but solve for the balance of ex- 
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probability meth od and the La r ge Velocity G r adient 
(LVG) method dSobolevI Il960t be Jong et al.1 1 1975 



Goldreich & Scovillelll976n 7 These 'intermediate-level 



methods require knowledge of molecular collisional 
data, whereas the previous 'basic-level' methods only 
required spectroscopic and dipole moment information. 
This extra requirement limits the use of these methods 
to some extent, because collisional data do not exist 
for all astrophysically relevant species. The advantage 
is that column density, kinetic temperature and volume 
density can be constrained, if accurate collision rates 
are known. As with rotation diagrams, this method can 
be used to compute synthet ic spectra to be compared 
with d ata with a x 2 statistic djansenlll995 ; Leurini et al. 
2004). 



The most advanced methods drop the local approxi- 
mation and solve for the intensities (or the radiative 
rates) as functions of depth into the cloud, as well as of 
velocity. Such methods are usually of the Accelerated 
Lambda Iteration (ALI) or Monte Carlo (MC) type, al- 
though hybrids also exist. The performance and con- 
ver gence of such programs ha ve recently been tested 
by 



Van Zadelhoff et al 



(2002). Using such programs 
one can constrain temperatur e, density, and velocity 
gradients within sources (e.g.. IVan der Tak et al1ll999t 
Tafalla et al.1 120021: Ijakob et alJ l2007h. and, ifenough 
observations are available, even molecular abundance 



profil e s (e.g., |V an der Tak et al.l l2000at ISchoier et al. 



2002: Maretetal. 2005). especially when coupled to 


chemical networks (e.g 


.. Dotv et al. 


2004; Evans et al. 


2005: Goicoechea et al. 


2006). 



This paper presents the public version of a radiative 
transfer code at the 'intermediate' level. The assump- 
tion of a homogeneous medium limits the number of 
free parameters and makes the program a useful tool 
in rapidly analyzing a large set of observational data, 
in order to provide constraints on physi cal conditions , 
such as density and kinetic temperature (lJansenl ll995). 
The program can be used for any molecule for which 
collisional rate coefficients are available. The input for- 
mat for spectroscopi c and collisional da ta is that of the 
LAMDA database (Scho ier et al.l l2005r^l where an on- 
line calculator for molecular line intensified, based on 
our program, can also be founcQ. 
The paper is set up as follows. Section [2] describes the 
radiative transfer formalism and introduces our notation 
of the key quantities. Section[3]describes the formalism 
which the program actually uses, and discusses its im- 
plementation. Section|4]compares the results of the pro- 
gram to those of other programs. The paper concludes 
in §[5] with suggested future directions of astrophysical 
radiative transfer modeling. 



2 http : //www . strw . leidenuniv . nl/~moldata 

3 Tn this naner. the, 'strength' of a line is an intrinsic ouantitv 



2. Radiative transfer and molecular 
excitation 

This section summarizes the formalism to analyze 
molecular line observations which our program adopts. 
For more detailed d isc ussions of radiative transfer see, 
e.g.. lCannonldl985l) or lRvbicki & Lightmanl dl979h . 



2.1. Basic formalism 

Describing the transfer of radiation requires a quantity 
which is conserved along its path as long as no local 
absorption or emission processes take place, and which 
includes the direction of travel. The quantity that sat- 
isfies this requirement is the specific intensity I v , de- 
fined as the amount of energy passing through a surface 
normal to the path, per unit time, surface, bandwidth 
(measured here in frequency units), and solid angle. The 
transfer equation for radiation propagating a distance ds 
can then be written as 



dl v 
ds 



Jv 



(1) 



where j v and a v are the local emission and extinction 
coefficients, respectively. The two terms on the right- 
hand side may be combined into the source function, 
defined by 

jv 



S v = 



(2) 



Writing the transport equation in its integral form and 
defining the optical depth, dr v = a v ds, measured along 
the ra)|f| one arrives at 



I v = I v (0)e- T *+ S v (T')e 



f 

Jo 



-(r r -r' y ) 



dr' (3) 



where l v is the radiation emerging from the medium and 
7y(0) is the 'background' radiation entering the medium. 
The above equations hold both for continuum radiation, 
which is emitted over a large bandwidth, and for spectral 
lines, which arise when the local emission and absorp- 
tion properties change drastically over a very small fre- 
quency interval, due to the presence of molecules. From 
this point the discussion will focus on bound-bound 
transitions within a multi-level molecule consisting of 
N levels with spontaneous downward rates A„i, Einstein 
coefficients for stimulated transitions B„i and Z?; u , and 
collisional rates C„i and C/„, between upper levels u and 
lower levels I. 

The rate of collision is equal to 



C„i = n co \y u i, 



(4) 



where n co \ is the number density of the collision partner 
(in ctrT 3 ) and y u \ is the downward collisional rate coeffi- 
cient (in cm 3 s _1 ). The rate coefficient is the Maxwellian 
average of the collision cross section, <x, 



Jul = 



-1/2 
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where E is the collision energy, k is the Boltzmann con- 
stant, 7\j n is the kinetic temperature, and fi is the re- 
duced mass of the system. The upward rates are ob- 
tained through detailed balance 



8i 



(6) 



where g,- is the statistical weight of level i. 

The local emission in transition u— >/ with laboratory 

frequency v„/, can be expressed as 



Jv 



hv u i 

47T 



(7) 



where «„ is the number density of molecules in level u 
and <py is the frequency-dependent line emission profile. 
The absorption coefficient reads 



hv u i 

a v = — — («; But <p Y 

4-7T 



■ n u B u iXv) , 



(8) 



where <p v and Xv are the line profiles for absorption and 
stimulated emission (counted as negative extinction), 
respectively. 

From here on we assume complete angular and fre- 
quency redistribution of the emitted photons, so that 
(p v =ip v =Xv, which is strictly only valid when collisional 
excitation dominates. This assumption allows the source 
function to be written as 



WkAh 



2hv li[g u ni 



niB lu - n u B ul 



— -1 , (9) 

gWu 



where we have used the Einstein relations. It is com- 
mon to introduce an excitation temperature T ex defined 
through the Boltzmann equation 

^ = exp [-(£„ - E,)/kT ex ] , (10) 
"/ gi 

where £, is the energy of level i, such that S Vlll = 
B v (T ex ), the specific intensity of a blackbody radiating 

at r ex . 

In the interstellar medium, the dominant line broadening 
mechanism is Doppler broadening. Except in very cold 
and dark cloud cores, observed line widths are much 
larger than expected from the kinetic temperature: this 
effect is commonly ascribed to random macroscopic gas 
motions or 'turbulence'. The result is a Gaussian line 
profile 



1 



V D V7T 



exp 



/ V «'\ 2 i 2 

-\v - v u , - v ■ n —J /v D 



(11) 



where vd is the Doppler width, v is the velocity vector 
of the moving gas at the position of the scattering, n is 
a unit vector in the direction of the propagating beam of 
radiation, and c is the speed of light. The Doppler width 
is the 1 /<? half-width of the profile, equal to A V/2 Vln2 
where AV is its full width at half-maximum. 
If the level populations n, are known, the radiative trans- 
fer equation can be solved exactly. In particular, under 



circumstellar media, the density is too low to attain LTE, 
but statistical equilibrium (SE) can often be assumed: 



dn 
dt 



i = o = n J p ji - Yj Pi ' =Ti ~ mDi 



(12) 



where the destruction rate coefficient of level i, and 
its formation rate coefficient Pj, are given by 



J Ay + BijJ v + Ci 

[Bjjjy + Cij 



a > d 
a < j). 



InEq.[l3l 



BijJ v 



B, 



Jo 



J Y (p( v) dv 



(13) 



(14) 



is the number of induced radiative (de-)excitations from 
state i to state j per second per particle in state i, and 



dO 



(15) 



is the specific intensity 7 V integrated over solid angle dQ. 
and averaged over all directions. The SE equations thus 
include the effects of non-local radiation. 
This discussion assumes that the state-specific rates of 
formation T t [cm 3 s~'] and destruction £), [s _1 ] are zero 
to ensure that the radiative transfer is solved indepen- 
dently of assumptions about chemical processes. In gen- 
eral, formation and destruction processes should be in- 
cluded explicitly to be able to deal with situations in 
which the chemical time scales are very short or the 
radiative lifetimes very long. For example, the forma- 
tion te mperature (in Tj) a ffects the rotational excitation 
of C 3 ( lRoueffetal.l[2 0Q2) and the vibrational excitation 
of H, Black & van pishoeckl ll987tl Burton et al.l ll990l: 
iTakahashi & Ueharall200ll) . systems for which line ra- 
diation only occurs as slow electric quadrupole transi- 
tions. The rotationa l excitation of reactive ions like CO + 
dFuente et al.ll2000t lBlacklll998h is also sensitive to Ti 
and T>i because the rates of reactions with H and H2 
rival the inelastic collision excitation rates. Similar con- 
siderations apply to the excitation of H3 in the Sgr A 
region close to the Galactic Center ( Van der Takl l2006). 
where electron impact excitation competes with disso- 
ciative recombination. 



2.2. Molecular line cooling 

Once the radiative transfer problem has been solved and 
the level populations are known, the cooling (or heating) 
from molecular line emission can be estimated. Since 
the level populations contain all the information of the 
radiative transfer, a general expression for the cooling is 
obtained from considering all possible collisional tran- 
sitions 

A = n ^ 2j X/"' r '" ~ n "7'»') /zy "'' ( 16 ) 

i l u>! 
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are in detailed balance at the kinetic temperature; there- 
fore it is possible for net heating to occur (A < 0) in 
cases where the crucial level populations have r ex >7\j n , 
owing to strong radiative excitation in a hot external ra- 
diation field. 

2.3. Escape probability 

The difficulty in solving radiative transfer problems is 
the interdependence of the molecular level populations 
and the local radiation field, requiring iterative solution 
methods. In particular, for inhomogeneous or geometri- 
cally complex objects, extensive calculations with many 
grid points are required. However, if only the global 
properties of an interstellar cloud are of interest, the cal- 
culation can be greatly simplified through the introduc- 
tion of a geometrically averaged escape probability ft, 
the probability that a photon will escape the medium 
from where it was created. This probability depends 
only on the optical depth t and is related to the inten- 
sity within the medium, ignoring background radiation 
and any local continuum, through 



(17) 



Several authors have developed detailed relations be- 
tween p and t for specific geometrical assumptions. 
Our program offers the user a choice of three such ex- 
pressions. The first is the expression derived for an ex- 
panding spherical shell, the so-called S obolev or large 
velocity grad i ent (LVG) app roximation dSobolevf l960: 
ICastorill97ri lElitzurl 1 1 992l p. 42-44). This method is 
also widely applied for moderate velocity gradients, to 
mi mic turbulent m otio ns. Our program uses the formula 
by lMihalasI d 1978b and lDe Jong et alJ dl980l) for this ge- 
ometry: 



1 

PlXG - - 



- T dr' = 



1 



(18) 



Second, in the case of a static, spherically symmet- 
ric and homogeneous medium the escape probability is 
dOsterbrock & Ferlandll2006l Appendix 2% 



/^sphe: 



L5 

T 



2/2 2 
- + - 



(19) 



Third, for a plane-parallel 'slab' geometry, applicable 
for instance to shocks, 



Alab = 



1 -, 



3t 



(20) 



is derived dDe Jong et al.lll975l) . Figure [T] plots the be- 
haviour of /? as a function of t f or these three cases; for 
more detail ed comparisons s ee IStutzki & Winnewisseil 
d 1985b and bssenkopf et al.1 d200lb . Users of our pro- 
gram can select either expression for their calculations. 
The on-line version of the program uses the formula for 
the uniform sphere, Eq. t!9) . 




Fig. 1. Escape probability j8 as a function of optical 
depth t for three different geometries: uniform sphere 
(solid line), expanding sphere (dotted line) and plane- 
parallel slab (dashed line). 

3. The program 

RADEX is a non-LTE radiative transfer code, written 
originally by J. H. Black, that uses the escape proba- 
bility formulation assuming an isothermal and homo- 
geneous medium without large-scale velocity fields. 
With the current increase of observational possibilities 
in mind, we have developed a version of this program 
which is suitable for public use. A guide for using the 
code in practice is provided in Appendix|A]and on-linfl 
Appendix [B] describes the adopted coding style. This 
section focuses on the implementation of the formalism 
of §|2]in the program. 



3.1. Basic capabilities 

For a homogeneous medium with no global velocity 
field, the optical depth at line centre can be expressed 
usingEqs. (f2l|7l|9l[IB, as 



A u iN mo \ 



8ttv 3 , 1.064AV 

lit 



gi 



(21) 



where N mo \ is the total column density, A V the full width 
at half-maximum of the line profile in velocity units, 
and xi the fractional population of level i. The formal- 
ism is analogous to the LVG method, with the global 
nj(dVjdK) re placed by the local N /AV, as in microtur- 
bulent codes dLeung & Liszt 1976). The program itera- 
tively solves the statistical equilibrium equations start- 
ing from optically thin statistical equilibrium (§ I3.41 i for 
the initial level populations. 

The program can handle up to seven collision partners 
simultaneously. In dense molecular clouds, H2 is the 
main collision partner for most species, but in some 
cases, separate cross sections may exist for collisions 
with the ortho and para forms of H2, and electron col- 
lisions may be important for ionic species. In diffuse 
molecular clouds and PDRs, excitation by atomic H be- 
comes important, particularly for fine structure lines, 
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refer to Flower d 19891) for the basi c theory of molecular 
collisions, and to iDubernetl ([2005) for an update of the 
latest results. 

The output of the program is the background-subtracted 
line intensity in units of the equivalent radiation temper- 
ature in the Rayleigh-Jeans limit. The background sub- 
traction follows traditional cm- and mm-wave spectro- 
scopic observations where the differences between on- 
source and off-source measurements are recorded, such 
that 



T R = 



2kv 2 



(/r-'v s ). 



(22) 



The radiation peak temperature Tr can be directly com- 
pared to the observed antenna temperature corrected 
for the optical efficiency of the telescope. However, it 
should be emphasized that RADEX contains no informa- 
tion about the geometry or length scale and that it is 
assumed that the source fills the antenna beam. If the 
source is expected to be smaller than the observational 
beam, computed line fluxes must be corrected before 
comparing to observed fluxes. 

In other types of observations, the continuum may 
not be subtracted from the data. In (sub-)millimeter 
and THz observations, for example with ESA's fu- 
ture Herschel space observatory, the dust continuum 
of many sources will be much stronger than any in- 
strumental error, and baseline subtraction may not be 
needed. The same is true for interferometer data, where 
the instrumental passband is well characterized. 

3.2. Background radiation field 

The average Galactic background (interstellar radia- 
tion field, ISRF) adopted in RADEX consists of several 
components. The main contribution is the cosmic mi- 
crowave background (CMB) whose absolute tempera- 
ture is taken to be Tcbr — 2.725+0.001 K based on the 



full C OBE data set as analyzed by iFixsen & M ather 
(2002). This model of the microwave background rep- 
resents the broadband continuum only and does not in- 
clude the strong emission lines, several of which contain 
significant power in the far-infrare d and (sub-)millime - 
ter part of the spectrum (see, e.g., IFixsen etaL I ll999l) . 
The ultraviolet/visible/near-infrared part of the spec- 
trum is based on the model of average Galactic sta rlight 
in the solar neig hborhood of MathisetalJ l[l983). The 
far-infrared and (sub-)millimeter part of the spectrum is 
based on the single-te mperature fit to the Galactic ther- 
mal dust emission o f Wri ght et al.l dl99ll) . At frequen- 
cies below 10 cm -1 (30 GHz), there is a background 
contribution from non-thermal radiation in the Galaxy. 
A tabulation of this spectrum in ASCII format is avail- 
able on- lindl and a graphical representation is shown in 
lBlackl(ll994 . 

One subtle aspect of the calculation is the distinction 
between the background seen by the observer and the 
background seen by the molecules. The continuum con- 
tribution to the rate equations may be composed of (1) 



region and (2) an internal continuum that arises within 
the emitting region. The CMB and ISRF are examples 
of external continuum components; dust emission from 
the line-emitting region is an example of an internal 
continuum. While an external continuum always fills 
the entire sky, an internal continuum may only fill a frac- 
tion of it, for example in the case of a circumstellar disk. 
With this distinction in mind, the internal intensity be- 
comes 

Jf = {3[B v (Tcbr) +T]I u v sa ] +(1 -P)lB v (T ex ) +0( 1-71)1^] 

(23) 

where /" sel is the continuous spectrum defined by the 
user. The factor rj, is the fraction of local continuum 
which arises outside the line emitting region, and the 
factor 9 is the fraction of local sky filled by the internal 
continuum. 

3.3. Chemical formation and destruction rates 

The equations of statistical equilibrium ( fT2l include 
source and sink terms. By default, RADEX sets the 
destruction rates equal to the same small value, 
Di = D = 10~ 15 s~', appropriate for cosmic -ray 
ionization plus cosmic-ray induced photod issociation 
dPrasad & Tarafdarll 19831; iGredel et al.lll989l) . The cor- 
responding formation rates are 

Ti = 10- 24 « tota ig,exp(-£,/^r fol - ra )/e(r fol - m ) (24) 

where n tota i is the sum of the densities of all collision 
partners, Tf 0rm is a formation temperature (default value 
300 K), and 

Q(T) = J] Si exp(-£ ; /£r form ) (25) 

;=I 

is the partition function. These assumptions imply a 
nominal fractional abundance of every molecule 



^mol 



= 10" 



(26) 



n total n totalf 

The value of the nominal abundance is inconsequential 
because the results in RADEX depend on N mo i/AV, but 
not on the fractional abundance. For most molecules 
currently in the associated database (LAMDA) and for 
the most commonly encountered interstellar conditions, 
these choices will not affect the observable excitation. 
The formation and destruction rates are computed in a 
subroutine that can be modified by the user to provide 
a more realistic description of chemical processes. For 
example, users may treat the combined ortho/para forms 
of molecules by introducing a realistic rf orm , especially 
in cases where no o/p interchange processes is likely 
to be effective. Other cases of potential interest include 
the photodissociation of large molecules into smaller 
molecules, or the evaporation of icy grain mantles into 
the gas phase. Our formulation in terms of a volume 
rate of formation is chosen to be independent of the de- 
tails of the formation process. In general, formation and 
destruction processes are important for molecules that 
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3.4. Calculation 

The input parameters of RADEX and its output are de- 
scribed in Appendix [A] Calculations with RADEX pro- 
ceed as follows. A first guess of the populations of the 
molecular energy levels is produced by solving statisti- 
cal equilibrium in the optically thin case. The only ra- 
diation taken into account is the unshielded background 
radiation field; internally produced radiation is not yet 
available. The solution for the level populations allows 
calculation of the optical depths of all the lines, which 
are then used to re-calculate the molecular excitation. 
The new calculation treats the background radiation in 
the same manner as the internally produced radiation. 
The program iteratively finds a consistent solution for 
the level populations and the radiation field. When the 
optical depths of the lines with t > 10 2 are stable from 
one iteration to the next to a given tolerance (default 
10~ 6 ), the program writes output and stops. 

3.5. Results 

There are several ways in which RADEX can be used 
to analyze molecular line observations. In most of 
these applications, the modeled quantity is the velocity- 
integrated line intensity, as the excitation is assumed 
to be independent of velocity. As a consequence, self- 
absorbed lines cannot be modeled satisfactorily with 
RADEX. In the simplest case, the temperature and den- 
sity are known from other observations and only the col- 
umn density of the molecule under consideration needs 
to be varied to get the best agreement with the observed 
line intensity. If the H2 column density is known from 
other observations, for example from an optically thin 
CO isotopic line, the ratio of the two column densi- 
ties gives the molecular abundance, averaged over the 
source. The RADEX distribution contains a Python script 
to automate this procedure which is further described in 
AppendixlDl 

Another often-used application of RADEX is to deter- 
mine temperatures and densities from the observed in- 
tensity ratios of lines of the same molecule. If the abun- 
dance of the molecule is constant throughout the source, 
the ratios should give source-averaged physical con- 
ditions independent of the specific chemistry of the 
molecule. Appendix ICl presents illustrative plots of line 
ratios for commonly observed molecules and lines in 
the optically thin case. For higher optical depths, the 
qualitative trends remain the same but there are quanti- 
tative differences. RADEX can readily be used to generate 
similar plots for moderately thick cases. Again, Python 
scripts are made available to automate this procedure 
(Appendix|D]). 

To illustrate the use of RADEX on actual observations, we 
take the observations of the HCO + 1-0 and 3-2 lines to- 
ward a relati vely simple source , the photon-dominated 
region IC 63 dJansen et all 1994 . The observed 1-0/3-2 
ratio corrected for beam dilution is 5.5+1.5. The kinetic 
temperature of the source is constrained from CO ob- 
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Fig. 2. Comparison of the predicted line strengths for 
the 10 lowest rotational transitions of CO for a homo- 
geneous isothermal sphere, with «e 2 = 10 5 cm~ 3 and 
7\in = 50 K, using different methods. Upper panel: The 
total optical depth through the sphere at line centre r 
and excitation temperature T sx as a function of the up- 
per rotational level J involved in the transition. Middle 
panel: The radiation temperature 7r obtained for each 
transition using RADEX with different prescriptions of 
the escape probability (3 and compare d with the resul t 
from the Monte-Carlo code (MCC) of lSchoierj (2000). 
Also shown are the results for optically thin emission in 
LTE. Lower panel: Tr obtained from RADEX compared 
with the results from MCC, 6T R = T™ cc - 7r. 



inferred column density from the absolute intensities is 
8 x 10 12 cm 2 , which, together with the overall H2 col- 
umn density of 5 xlO 21 ctrT 2 , gives anHCO + abundance 
with respect to H 2 of 1.6 x 10~ 9 . 

A slightly more c omplicated situation aris es for the 
Orion Bar PDR dHogerheijde et all 1 19951) . For this 
source, both HCO + 1-0, 3-2 and 4-3 lines have been 
observed. The 1-0/3-2 ratio gives an order of magni- 
tude lower density than the 3-2/4-3 ratio. This differ- 
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density tracers such as CO are usually much less sensi- 
tive to density. One possible solution is a clumpy PDR 
model in which the 1-0 line is mostly produced in the 
low-density interclump gas containing 90% of the mate- 
rial and the 4-3 line in the high-density clumps. Within 
this clumpy model, a single column density fits all three 
lines and an accurate abundance can be derived. Note 
that this technique of adding the results of two models 
is only applicable at low optical depths. 
If only two lines of a molecule have been observed, 
the line ratio can be used as indicator of tempera- 
ture or density, depending on molecule and transition 
(Appendix [C]>- A single line ratio is never enough 
to constrain both temperature and density, though. 
For multi-line observations, a comparison of data and 
models in terms o f y 2 is preferred. See for ex ample 



Van d er Tak etail ( 2000bl). ISchoier et al] |2002), and 



Leurini et a l. (2004) for details of such calculations. 



3.6. Limitations of the program 

The current version of the program does not in- 
clude a contribution from continuous (dust or free-free) 
opacity to the escape probability, as for example in 
Takahashi et al.1 (119831) . Continuum radiation from dust 
is generally negligible at long wavelengths ( ! 1 mm) 
but becomes important for regions with very high col- 
umn densities (such as protoplanetary disks) and at 
far-infrared and shorter wavelengths (Z lOOyum). Free- 
free radiation may become important for the calcula- 
tion of atomic fine structure lines from H II regions; 



other programs such as CLOUDY (Ferland 20031J may 



be more suitable for this purpose. The absence of con- 
tinuous opacity limits the applicability of the program 
particularly in situations where infrared pumping is im- 
portant, either directly throug h rotational transitions or 
via vibrational transiti ons ( Carroll & Goldsmith! 1981 
Hauschildt et"aDl993l) . 



Another limitation of the program is that only one 
molecule is treated at a time, so that the effects of line 
overlap are not taken into account. Such overlaps may 
occur both at radio and at infrared wavelengths (e.g., 
Exposito et al . 2006). In special cases, overlap between 
lines of the same molecule may influence their excita- 
tion, f or example the hyp erfine components of HCN or 
N 2 H + JPaniel et al.ll2006l) . 

For certain molecules under certain physical conditions 
(especially low density and/or strong radiation field), 
population inversions occur, which cause negative op- 
tical depth and hen ce nonlinear amplification of the 
incoming radiation (Elit zurlll992l) . This phenomenon, 
known as 'maser' action, requires non-local treatment 
of the radiative transfer, in particular a fine sampling of 
directions, for which RAD EX is not set up. Generally, the 
escape probability approximation is justified until the 
masers saturate, which occurs at t m -1. In practice, 
the computed intensities of lines with tZ - 0.1 are not 
as accurate as those of other lines, and the intensities of 
lines with t < — 1 should be disregarded altogether Tf 
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Fig. 3. Excitation temperature of the HCO + 1-0 tran- 
sition as a function of radius for the model cloud of 
§ 14.1.21 calculated with RADEX assuming static spher- 
ical, expanding spherical, and slab geometry (dashed 
I dash-dotted / dotted lines), with a multi-zone es- 
cape probability program (long/short dashes) and with 
a Monte Carlo code (solid lines). The panels are for dif- 
ferent column densities, hence optical depths. Note the 
different vertical scales. 



of non-maser lines may also be affected. While special- 
ized programs should be used to calculate the intensities 
of maser lines (e.g.. ISp aans & v an Langeveldd 1 1 992 : 



Gray & Field! 19951: lYates et al.ll 19971) . RADEX may well 
be used to predict which lines of a molecule may display 
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which are pumped by infrared radiation (ILeurini et al 
2004 . 



4. Comparison with other methods 

This section shows a comparison of RADEX with other 
programs, first for the case of constant physical con- 
ditions (§ 14. Il l and second for variable conditions 
(§ 14.2b . Comparison is with the analytical rotation di- 
agram method and with Monte Carlo methods, which 
have been benchm arked to high accuracy , both for 
the ca se of HCO + jVan Zadelhoff et alJ kOOaH and of 



H 2 (I Van der Tak et all 1200511 u l Throughout this sec- 
tion, mole cular data have bee n taken from the LAMDA 
database (ISchoier et al ] l2005h . 



4.1. Homogeneous models 

4.1.1. The case of CO 

To test the RADEX code, we have compared its output 
both to an optically thin LTE analysis (rotation diagram 
method) and a full r adiative transfe r analysis using a 
Monte-Carlo method (ISchoier 2000). The test problem 
consists of a spherically symmetric cloud with a con- 
stant density, nCEL;), of 1 x 10 5 ctrT 3 within a radius of 
100 AU. In this example only the CO emission is treated 
using a fractional abundance of 1 x 10~ 4 relative to H2 
yielding a central CO column density of Nco = 3 x 
10 16 cirT 2 and an average value of Af=2xl0 16 ctrT 2 . The 
kinetic temperature is set to 50 K, the background tem- 
perature to 2.73 K, and the line width to AV= 1 .0 km s . 
Fig. [2] presents the results of the calculations for the 
ten lowest rotational transitions. The excitation temper- 
atures of the lines vary from being close to thermalized 
for transitions involving low /-levels, to sub-thermally 
excited for the higher-lying lines. The optical depth in 
the lines is moderate (~ 1 - 2) to low. It is seen that the 
expressions of the escape probability for the uniform 
sphere and the expanding sphere give almost identical 
solutions which are close to that obtained from the full 
radiative transfer (MCC in Fig. The slab geometry 
gives slightly higher intensities, in particular for high- 
lying lines. The optically thin approximation, where the 
gas is assumed to be in LTE at 50 K, produces much 
larger discrepancies, up to a factor of ~ 2, and only gives 
the correct intensity for the J = 1 — > line, where the 
LTE conditions are met. 

4.1.2. The case of HCO+ 

To further verify the performance of the RADEX 
program, we have compared its results to that 
of another program that does not use the local 
approximation: the Monte Carlo program RATRAN 
(iHogerheijde & van der Takl |2000). We also compare 
the results to thos e from the multi-zone esc ape proba- 
bility program by iPoelman & Spaansl d2005l) . The test 



case is a cloud with n(H2) = lxl0 4 cirT 3 , 7kj n =10K, 
7b g =2.73K, and a line width of A V= 1.0 km s _1 , equiv- 
alent to b/)=0.6km s _1 . The pure rotational emission 
spectrum of HCO + was calculated for column densi- 
ties of 10 12 , 10 13 , 10 14 and 10 15 cirT 2 , which for RADEX 
were given directly as input parameters. For the multi- 
zone programs, a cloud radius of 10 18 cm was specified 
along with abundances of 10~ 10 -10~ 7 , distributed over 
50 cells. 

Figure [3] shows the calculated excitation temperature 
of the HCO + 1-0 transition as a function of radius 
for these physical conditions. For A^(HCO + ) 1 10 12 cm" 2 , 
the excitation is independent of radius and the calcu- 
lations for the various geometries agree to »T0%. The 
dependence of the excitation on radius and on geom- 
etry increases with increasing column density, and for 
N(HCO + );10 15 cirr 2 , the curvature of the T ex distri- 
bution becomes too large to ignore. The corresponding 
line optical depth is »T00, with ^20% spread between 
the various estimates (Fig.HJl. The curvature arises be- 
cause at the cloud center, photon trapping thermalizes 
the excitation, while at the e dge, the emission can es- 
cape the cloud (Bernesfl979l) . We do not recommend to 
use RADEX at line optical depths ! 100, because the cal- 
culated excitation temperature may not be representa- 
tive of the emitting region. However, even if some lines 
are highly optically thick, RADEX may well be used to 
analyze other lines which are optically thin. For exam- 
ple for H2O, the ground state lines often have t ~ 1000, 
but RADEX is well capable of computing intensities for 
higher-lying transitions which are not as optically thick. 
At low optical depth, variations in T ex translate directly 
into changes in emergent line intensity. Thus, differ- 
ences as large as 20% in calculated line flux can arise 
depending on the choice of escape probability descrip- 
tion, even for moderately thick cases. At high optical 
depth, the direct connection between T ex and line flux 
is lost because of the dependence on the adopted veloc- 
ity field. The assumption in the program that the opti- 
cal depth is independent of velocity breaks down in this 
case. In this limit, the peak line temperature Tr gives 
the value of T ex at the t = 1 surface of the cloud in this 
specific transition. 

The results shown in this section do not translate easily 
to other HCO + lines such as 7=3— >2, because the ex- 
citation is governed by several competing effects. The 
optical depth of the J-3— >2 line may be higher or lower 
than that of the /=1— »0 line, depending on temperature 
and density. Observers are encouraged to use RADEX to 
study the excitation of their lines as a function of these 
parameters, and also consider geometric variations. 



4.2. Observations of a Young Stellar Object 

To compare a typical RADEX analysis with other meth- 
ods for a situation which varying physical conditions, 
we choose a molecule for which many lines can be 
observed: para-formaldehyde, P-H2CO. Figure[5]shows 
observations of n-HoCO fsnb-imillimeter emission lines 
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Fig. 4. Optical depth of the HCO + 1-0 transition for 
the model of § 14.1.21 calculated with RADEX assuming 
static spherical, expanding spherical, and slab geometry 
(filled triangles/filled squares/open triangles), with the 
multi-zone escape probability code (open circles) and 
with the Monte Carlo code (filled circles). 
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Fig. 5. Line strengths of P-H2CO observed towards 
the embedded low-mass protostar IRAS 16293-2422 
(squares with error bars), modeled assuming LTE (solid 
line), using RADEX (triangles), and using a Monte Carlo 
program assuming a constant abundance (solid circles) 
and an abundance varying with radius (open circles). 
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Fig. 6. Distributions of the x 2 parameter correspond 
ing to the models in Figure [3] Th e RADEX results are 
for »(H 2) = 10 6 cirT 3 as found by IVan Dishoeck et al 
(1995). 



Van Dishoeck et al. (1995). The data are analyzed us- 



ing three methods: assuming LTE (with a rotation dia- 
gram), assuming SE (using RADEX), and using a Monte 
Carlo program. The free parameters for the LTE fit are 
the excitation temperature T ex and the column density 



Figure[6]shows the distributions of the^ 2 parameter for 
the LT E and SE fits, calculated in the standard way (see, 
e.g.. IVan der Tak et alJEoOOri ISchoier et al.ll2002l) as- 
suming a 20% uncertainty for all observed points except 
the line with the highest upper level energy where a 30% 
uncertainty was used. As seen from the figure, the non- 
LTE method gives a better fit to the data, as quantified 
by the lower minimum y 1 value. This result is not nec- 
essarily surprising given that more free parameters are 
available. A more important difference is that the esti- 
mates of temperature and column density between the 
two methods are substantially different, in particular the 
temperature (50 vs 150 K). Since the non-LTE method 
involves fewer assumptions about the physical state of 
the cloud, its results are to be preferred. 
These results illustrate that rotation diagrams may give 
misleading results when determi ning physical proper- 



ties of interstellar gas clouds (cf. Ijohnstone et al. 2003 



for the case of CH3OH). Figure|5]also demonstrates that 
temperatures and column densities derived from rota- 
tion diagrams tend to depend on which lines happen to 
have been observed (cf. the HCO + case in § l4.1.2t . From 
other data, IRAS 16293-2422 is actually known to have 
a gradient in temperature and density throughout its en- 
velope, which cannot be modelled properly with either 
technique. For such situation, a full Monte Carlo radia- 
tive transfer method is needed in which both the physi- 
cal conditions and the abundances can vary with radius 
(Fig- El circles). Nevertheless, the column densities and 
abundances inferred with RADEX using the physical con- 
ditions inferred from the line ratios differ by only a fac- 
tor of a few from those found with the more sophis- 
ticated analysis, at least for the particula r zone of the 



sourc e to which those conditions apply dSchoier et al 



2002). 



5. Conclusions 

We have presented a computer program to analyze spec- 
tral line observations at radio and infrared wavelengths, 
based on the escape probability approximation. The pro- 
gram can be used for any molecule for which collisional 
data exist; such input data are available in the required 
format from the LAMDA database. The program can be 
used for optical depths from «-0. 1 to a;100. 
The limited number of free parameters makes RADEX 
very useful to rapidly analyze large datasets. As an 
example, observed line intensity ratios may be com- 
pared with the plots in Appendix ICl to estimate density 
and kinetic temperature. Ratios of other lines and other 
molecules may be easily computed using the Python 
scripts included in the RADEX distribution. The program 
may also be used to create synthetic spectra. This capa- 
bility will be important to model the THz line surveys 
from the HIFI instrument onboard the Herschel space 
observatory. 

In the future, we plan to inc orporate a multi-zone es- 



cape probability formalism dPoelman & Spa ans 2005; 
Elitzur & Asensio Ramosl I2QQ60 which will enable 
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depths, the calculation may also start from LTE condi- 
tions rather than from optically thin statistical equilib- 
rium. Robust convergence may be achieved by starting 
from either initial condition and requiring the two an- 
swers to be equal. For the modeling of crowded spec- 
tra, the effects of line overlap will also need to be con- 
sidered, for instance in the 'all or nothing' approach 
(Cesaroni & Walmslevl[l99lh . Such spectra will be rou- 
tinely observed with the superb resolution and sensitiv- 
ity of ALMA. 

Our program is free for anybody to use for science, pro- 
vided that appropriate reference is made to this paper. 
For any other purpose such as to incorporate the pro- 
gram into other packages which may be distributed to 
the public, prior agreement with the authors is needed. 
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Appendix A: Program input and output 

A.1. Program input 

The input parameters to RADEX are the following: 

1 . The name of the molecular data file to be used. 

2. The name of the file to write the output to. 

3. The frequency range for the output file [GHz]. All 
transitions from the molecular data file are always 
taken into account in the calculation, but often it is 
practical to write only a limited set of lines to the 
output. 

4. The kinetic temperature of the cloud [K]. 

5. The number of collision partners to be used. Most 
users will want H2 as only collision partner, but in 
more specialized cases, additional collisions with 
H or electrons may for instance play a role. See 
the molecular datafiles for details. For some species 
(CO, atoms) separate collision data for ortho and 
para H2 exist; the program then uses the thermal or- 
tho/para ratio unless the user specifies otherwise. 

6. The name (case-insensitive) and the density [cm" 3 ] 
of each collision partner. Possibilities are H2, P-H2, 
0-H2, electrons, atomic H, He, and H + . 

7. The temperature of the background radiation field 
[K]. 

- If >0, a black body at this temperature is used. 
Most users will adopt the cosmic microwave 
background at rcMB=2.725(l+z)K for a galaxy 
at redshift z. 

- If =0, the average interste llar radiation field 
(ISRF) is used, taken from iBlackl (1 1994 with 
modifications described in § 13.21 This spectrum 
is not adjustable by a scale factor because it con- 
sists of several components that are not expected 
to scale linearly with respect to each other. 

- If <0, a user-defined radiation field is used, spec- 
ified by values of frequency [cm -1 ], intensity 
[Jy nsr~'], and dilution factor [dimensionless]. 
Spline interpolation and extrapolation are ap- 
plied to this table. The intensity need not be 
specified at all frequencies of the line list, but 
a warning message will appear if extrapolation 
(rather than interpolation) is required. 

8. The column density of the molecule [cm -2 ]. 

9. The FWHM line width [km s -1 ]. 

A.2. Program output 

The output file written by RADEX first replicates the in- 
put parameters, and then lists the following quantities 
for each spectral line within the specified frequency 
range. 

1. Quantum numbers, upper state energy [K], fre- 
quency [GHz], and wavelength [/mi]. These num- 
bers are just copied from the molecular data file, 
which usually comes from the LAMDA database. 
Frequencies from this database are generally of 



line catalogs such as CDMS dMiiller et al.l l200lF 
should be consulted. 

2. The excitation temperature [K] as defined in 
Eq. ( [Tol l. In general, different lines have different 
excitation temperatures. Lines are thermalized if 
T e x-Tyin, in LTE, all lines are thermalized. 

3. The line optical depth, defined as the optical depth 
of the equivalent rectangular line shape (<p v = 1/Av). 

4. The line intensity, defined as the Rayleigh-Jeans 
equivalent temperature Tr [K] . 

5. The line flux, defined as the velocity-integrated 
intensity, both in units of Kkm s _1 (common in 
radio astronomy) and of ergcirT 2 s~' (common 
in infrared astronomy). The line flux is calcu- 
lated as 1.0645r*AV, where the factor 1.0645 = 
sfn/ (2 Vln 2) converts the adopted rectangular line 
profile into a Gaussian profile with an FWHM of 
AV. The integrated profile is useful to estimate the 
total emission in the line, but it has limited mean- 
ing at high optical depths, because the change of 
optical depth over the line profile is not taken into 
account. Proper modeling of optically thick lines re- 
quires programs that resolve the source both spec- 
trally and spatially (see § 14.1.21 for further discus- 
sion). 

Auxiliary output files can be generated, for example to 
display the adopted continuum spectrum. 

Appendix B: Coding standards 

The original version of RADEX was written in such a 
way as to minimize the use of machine memory which 
was expensive until a decade ago. Nowadays, clarity 
and easy maintenance are more important requirements, 
which is why the source code has been re-written fol- 
lowing the rules below. We hope that these rules will be 
useful for the development of other 'open source' as- 
tronomical software. For further guidelines on scientific 
programming we recommend the Software Carpentri* 3 ! 
on-line course. 

1. All the action is in subroutines; the sole purpose of 
the main program is to show the structure of the pro- 
gram. The subroutines are grouped into several files 
for a better overview; compilation instructions for 
automated builds on a variety of platforms are in a 
Makefile. 

2. The program text is interspersed with comments at 
a ratio of * 1:1. In particular, each subroutine starts 
with a description of its contents, its input and out- 
put, where incoming calls come from, and which 
calls go out. Then the properties of each variable are 
described: contents, units and type. 

3. Variables and subroutines have descriptive names 
with a length of 5-10 letters. Names of integer vari- 
ables start with the letters i..n; names of floating- 
point variables (always of double precision) with 
a..h or o..z. There are no specific namings for vari- 
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in CLOUDY. Names are always based on the English 
language. We do not use upper-, lower-, or camel- 
case to distinguish types of identifiers, as some pro- 
grams do. 

4. Loops are marked by indenting the program text. 
The loop variables are always called ilev, iline .., 
never just i. 

5. Subroutines start with a check whether the input pa- 
rameters have reasonable values. Such checks force 
soft landings if necessary, and avoid runtime errors. 

6. Statements with calculations use spaces around 
the =, + and - symbols, but not around others. 
Calculations that consist of multiple steps are split 
over as many program lines. Multiple assignments 
in a row are aligned at the = sign. 

7. The program text avoids "magic numbers" both in 
calculations and in definitions. Numbers that are of- 
ten used such as Planck's constant are defined at one 
central place in the program. Similarly, often-used 
variables such as physical parameters are stored in 
shared memory rather than passed on via subroutine 
calls. 



Appendix C: Diagnostic plots of molecular 
line ratios 

We have used the RADEX program and the LAMDA 
database to calculate line ratios of several commonly 
observed molecules for a range of kinetic tempera- 
tures and H2 densities. The plots in this Appendix may 
be used by observers to estimate physical conditions 
from their data. Line ratios have the advantage of be- 
ing less sensitive to calibration errors than absolute line 
strengths, especially when the two lines have been mea- 
sured with the same telescope, receiver and spectrome- 
ter. 

The calculations assume a column density of 10 12 cm -2 , 
a line width of 1.0 km s _1 , and use a 2.73 K blackbody 
as background radiation field. Under these conditions, 
the lines are optically thin, so that the line ratios do not 
depend on column density. The calculations also assume 
that the emission in both lines fills the telescope beams 
equally, which may be the case if the lines are close 
in frequency. However, lines are generally measured in 
beams of different sizes, and the observations need to 
be corrected to account for this effect, if the source is 
known to be compact. 

Linear molecules such as CO (Fig. IC 1 b are tracers 
of density at low densities, when collisions compete 
with radiative decay. At higher densities, the excitation 
becomes thermalized and the line ratios are sensitive 
to temperature. For a given molecule, moving up the 
/-ladder means probing higher temperatures and densi- 
ties. Note that for the column densities of typical dense 
interstellar clouds, the CO lines are optically thick, and 
observations of 13 CO or even rarer isotopologues must 
be used to probe physical conditions. 
The critical densities of molecular lines scale as p^v 3 , 
where p is the permanent dipole moment of the 
molecule and v is the frequency of the line. Indeed, the 
CS molecule (~Fig. IC.2b has a larger dipole moment than 
CO, and its line ratios are mainly probes of the den- 
sity. The small frequency spacing between the lines of 
CS makes this molecule very useful to probe density 
structure (e.g., IVan der Tak et all l2000bl) . The HCO + 
and HCN molecules (Fig. IC.3t display similar trends to 
CS, although their line spacing is not as small. 
Non-linear molecules such as H2CO (Figs. lC~4l and lC.5l l 
have the advantage that both temperature and density 
may be probed within the same frequency range. Ratios 
of lines from different /-states tend to be density tracers 
(left panels), while ratios of lines from the same /-state 
but different /f-states are mostly temperature probes 
(right panels). The lines of H2CO are often quite strong, 
making this molecule a favourite tracer of temperature 
and density dMangum & Woottenlll993l) . Other asym- 
metric molecules have also bee n used, such as H2 CS 
dBlakeet all 1994b and CH 3 OH dLeurini et alj|2004b al- 
though abundance variations from source to source or 
even within sources of ten complicate the interpretation 
dJohnstone et alJl2003b . 
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Appendix D: The Python scripts 

The RADEX distribution comes with two scripts, 
radex_line .py and radex_grid.py, to automate 
standard modeling procedures. The scripts are written 
in Python and are run from the Unix shell command 
line after manual editing of parameters. 

The first script, radex_line .py, calculates the column 
density of a molecule from an observed line intensity, 
given estimates of kinetic temperature and H2 volume 
density. The input parameters are: 

1 . The kinetic temperature [K] 

2. The number density of H2 molecules [cm 4 ] 

3. The temperature of the background radiation field 
[K], usually 2.73 (CMB). 

4. The name of the molecule (or molecular data file) 

5. The frequency of the line [GHz] 

6. The observed line intensity [K] 

7. The observed line width [km s _1 ] 

Furthermore, two numerical parameters have good de- 
fault values but will need to be changed occasionally: 

1. The free spectral range around the line (default 
10%): this number must be smaller for molecules 
with many line close in frequency, such as CH3OH. 
The program uses this parameter to find the ob- 
served line from the list of lines in the molecular 
model. 

2. The required accuracy (default 10%): The default 
corresponds to the calibration uncertainty of most 
telescopes. 

The script iterates on column density until the observed 
and modeled line fluxes agree to within the desired ac- 
curacy. The best-fit column density is directly written 
to the screen. The file radex . out gives details of the 
best-fit model. 

The second script, radex_grid.py, runs a series of 
RADEX models to estimate the kinetic temperature 
and/or the volume density from an observed line ratio. 
The user needs to set the following input parameters: 

1. The grid boundaries: minimum and maximum ki- 
netic temperature [K] and minimum and maximum 
H2 volume density [cm -3 ]. 

2. The temperature of the background radiation field 
[K], usually 2.73 (CMB). 

3. The molecular column density [cm -2 ]. For the 
illustrative plots in Appendix [C] a low value 
(N-IQ 12 cirT 2 ) was used, so that the line ratios are 
independent of column density (optically thin limit). 
However, in modeling specific observations, it is 
worth varying this parameter to assess the sensitivity 
of the line ratio to column density. 

4. The observed line width [km s _1 ], usually an aver- 
age of the widths of the two lines. 



1 . The number of grid points along the temperature and 
density axes. 

2. The free spectral range around the line (see above) 

The name of the molecule and a list of observed line 
ratios and names of associated output files are given at 
the start of the main program. The script produces a file 
radex. out which is a tabular listing of temperature, 
log density, and line ratio. This results may be plotted 
with the user's favourite plotting program. 
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Fig. C.l. Line ratios of CO in the optically thin limit 
as a function of kinetic temperature and H2 density. 
Contours are spaced linearly and some contours are la- 
beled for easy identification. 
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Fig. C.2. Line ratios of CS in the optically thin limit 
as a function of kinetic temperature and H2 density. 
Contours are spaced linearly and some contours are la- 
beled for easy identification. 
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Fig. C.3. Line ratios of HCO + and HCN in the optically 
thin limit as a function of kinetic temperature and H2 
density. Contours are spaced linearly and some contours 
are labeled for easy identification. 
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Fig. C.4. Line ratios of 0-H2CO in the optically thin 
limit as a function of kinetic temperature and H2 den- 
sity. Contours are spaced linearly and some contours are 
labeled for easy identification. 
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Fig. C.5. Line ratios of P-H2CO in the optically thin 
limit as a function of kinetic temperature and H2 den- 
sity. Contours are spaced linearly and some contours are 
labeled for easy identification. 



